##################################################################
# This file runs the main and supplementary analyses.
##################################################################
#set working directory
setwd("/Users/agunderson/Dropbox/Gender and Policing")
#clear working environment
rm(list=ls())
#load packages
library(plyr)
library(dplyr) 
library(ggplot2)
library(tidyverse)
library(lfe)
library(stargazer)
library(xtable)
library(readxl)
library(cdlTools)
library(MatchIt)
library(cobalt)
options(scipen=999) #to allow merge in of very large numbers

##################################################################
# Model function
##################################################################
main_model <- function(dataframe,tabletitle,label,dv1,dv2,iv1,ivlab,depvarlabel){
  model1 <- felm(dv1~iv1
                 |year|0|0,data=dataframe)
  
  #model 2: female police officers and controls
  model1controls <- felm(dv1~iv1+totalemployees+log(budget)+percblackfulltime+raperate+sheriff+communitypolicing
                         |year|0|0,data=dataframe) 
  
  #Model 2: female police leadership with diff controls 
  model2controls <- felm(dv1~iv1+totalemployees+ percblackfulltime+raperate+sheriff
                         |year|0|0,data=dataframe)
  #with reports
  model1a <- felm(dv2~iv1|year|0|0,data=dataframe)
  
  #model 2: female police officers and controls
  model1acontrols <- felm(dv2~iv1+totalemployees+log(budget)+percblackfulltime+sheriff+communitypolicing
                          |year|0|0,data=dataframe) 
  
  #Model 2: female police leadership with diff controls 
  model2acontrols <- felm(dv2~iv1+totalemployees+percblackfulltime+sheriff|year|0|0,data=dataframe)
  
  #Export for Latex -----
  stargazer(model1a,model1,model2acontrols,model2controls,model1acontrols,model1controls,covariate.labels = c(ivlab,
                                                                                                              "Total Agency Employees","Logged Agency Budget",
                                                                                                              "Perc. Black Officers","Rape Report Rate","Sheriff",
                                                                                                              "Community Policing"),
            label=label,
            title=tabletitle,
            dep.var.labels=depvarlabel,
            style="apsr",notes="Year dummy variables included.",
            out= paste0("./Female Cops Paper/Tables/",label,".tex"),
            font.size = "tiny",table.placement="h!",column.sep.width = "-10pt")
}

year_fun_arrests_pre <- function(data){
  model_arrests <- felm(sumrapearrestsrate~percfemale+
                          totalemployees+
                          percblackfulltime+raperate+sheriff
                        |0|0|0,data=data) 
} 

year_fun_arrests_post <- function(data){
  model_arrests <- felm(sumrapearrestsrate~percfemale+
                          totalemployees+log(budget)+
                          percblackfulltime+raperate+sheriff+communitypolicing
                        |0|0|0,data=data) 
} 

year_fun_arrests_post2007 <- function(data){
  model_arrests <- felm(sumrapearrestsrate~percfemale+
                          totalemployees+log(budget)+
                          percblackfulltime+raperate+sheriff+communitypolicing+dvunit
                        |0|0|0,data=data) 
} 

year_fun_reports_pre <- function(data){
  model_reports <- felm(raperate~percfemale+
                          totalemployees+
                          percblackfulltime+sheriff
                        |0|0|0,data=data)
} 

year_fun_reports_post <- function(data){
  model_reports <- felm(raperate~percfemale+
                          totalemployees+log(budget)+
                          percblackfulltime+sheriff+communitypolicing
                        |0|0|0,data=data)
} 

year_fun_reports_post2007 <- function(data){
  model_reports <- felm(raperate~percfemale+
                          totalemployees+log(budget)+
                          percblackfulltime+sheriff+communitypolicing+dvunit
                        |0|0|0,data=data)
} 

##################################################################
# Import data --------------
################################################################## 
genderpolice <- read_csv("./Female Cops Paper/PoP/PoPData.csv")%>%
  mutate(rapeclearancerate = (rapescleared*100000)/population) %>%
  na_if(., Inf) %>%
  group_by(ORI) %>%
  filter(n()<10) %>% #some mismatches in names, so need to exclude those from analysis (more than one agency matched on name only)
  filter(percblackfulltime<101)%>% #two agencies have more black employees than total employees
  as.data.frame() 

genderpolice <- genderpolice %>% filter(sumrapearrestsrate<1000) #filter out massive outlier

##################################################################
# Figures --------------
################################################################## 
#figure 1
genderpolicenona <- genderpolice %>%
  dplyr::select(year,percfemale,totalemployees,percblackfulltime,raperate,sheriff,sumrapearrestsrate) %>%
  filter(!is.na(raperate)) %>%
  filter(!is.na(percblackfulltime)) %>%
  filter(!is.na(sumrapearrestsrate)) %>%
  mutate(largedepartment = if_else(totalemployees>155,1,0),
         highraperate = if_else(raperate>25,1,0),
         highfemale = if_else(percfemale>8,1,0))  

mod_match <- matchit(highfemale ~ totalemployees+sheriff+year, data = genderpolicenona,
                     method = "nearest")  

love.plot(mod_match,binary="std",var.names = data.frame(old = c("distance", "totalemployees", "sheriff","year"),
                                                        new = c("Propensity Score", "Total Employees", 
                                                                "Sheriff","Year")),colors=c("grey","black"))
ggsave("./Female Cops Paper/Figures/matchedplot.pdf") 

#figure 2
#see here for replication code: https://codeocean.com/capsule/8907164/tree/v1
#we ran as-is, except changed line 256 in s00_preprocess to:
#arrests <- fread(file.path('arrests.csv.gz')) %>% filter(., grepl('SEX|PORNOGRAPHY|OBSCENE|OBSCENITY|INDECENCY|INDECENT|PEEPING TOM|EXPOSE|CHILD LURING|LEWD|VIDEOTAPE THROUGH CLOTHES', statute_description))
#to only include arrests for sex crimes following CPD's guidelines (https://gis.chicagopolice.org/pages/crime_details)

#################################################################
# Tables
##################################################################
#table 1
genderpolice.cor <- genderpolice %>%
  ungroup()%>%
  select_if(is.numeric)%>%
  dplyr::select(sumrapearrestsrate,percfemale,totalemployees,percblackfulltime,raperate,sheriff,budget,communitypolicing)

cor <- as.data.frame(cor(genderpolice.cor,use="p"))

colnames(cor) <- c("Rape Arrest Rate","% Female Sworn Officers",
                   "Total # Employees",
                   "% Black Officers","Rape Report Rate","Sheriff","Budget","Community Policing")
rownames(cor) <- c("Rape Arrest Rate","% Female Sworn Officers",
                   "Total # Employees",
                   "% Black Officers","Rape Report Rate","Sheriff","Budget","Community Policing")

print(xtable(cor,label = "correlation",caption="Correlation Matrix"),
      file="./Female Cops Paper/Tables/correlation.tex", size="\\tiny")

#table 2
main_model(genderpolice,"The Effect of Female Police on Rape Report and Arrest Rates, 1987-2016: Using the Number of Female Officers",
           "femalecopsnumber",genderpolice$sumrapearrestsrate,genderpolice$raperate,genderpolice$femalefulltime,
           "No. Female Officers",
           c("Rape Report Rate","Rape Arrest Rate","Rape Report Rate",
             "Rape Arrest Rate","Rape Report Rate","Rape Arrest Rate"))

#table 3
model1 <- felm(sumrapearrests~percfemale
               |year|0|0,data=genderpolice)
model1controls <- felm(sumrapearrests~percfemale+
                         totalemployees+log(budget)+
                         percblackfulltime+raperate+sheriff+communitypolicing
                       |year|0|0,data=genderpolice) 
model2controls <- felm(sumrapearrests~percfemale+
                         totalemployees+
                         percblackfulltime+raperate+sheriff
                       |year|0|0,data=genderpolice)
stargazer(model1,model2controls,model1controls,covariate.labels = c("Perc. Female Officers",
                                                                    "Total Agency Employees","Logged Agency Budget",
                                                                    "Perc. Black Officers","Rape Report Rate","Sheriff",
                                                                    "Community Policing Component"),
          label="femalecopsnumbersexarrests",
          title="The Effect of Female Police on Rape Arrests, 1987-2016",
          dep.var.labels="Number of Rape Arrests",
          style="apsr",notes="Year dummy variables included.",out="./Female Cops Paper/Tables/femalecopsnumbersexarrests.tex",
          font.size = "small",table.placement="h!")

#table 4
model1 <- felm(sumarrestsrate~percfemale
               |year|0|0,data=genderpolice)
model1controls <- felm(sumarrestsrate~percfemale+
                         totalemployees+log(budget)+
                         percblackfulltime+raperate+sheriff+communitypolicing+violentcrimerate+propertycrimerate
                       |year|0|0,data=genderpolice) 
model2controls <- felm(sumarrestsrate~percfemale+
                         totalemployees+
                         percblackfulltime+raperate+sheriff+violentcrimerate+propertycrimerate
                       |year|0|0,data=genderpolice) 
stargazer(model1,model2controls,model1controls,covariate.labels = c("Perc. Female Officers",
                                                                    "Total Agency Employees","Logged Agency Budget",
                                                                    "Perc. Black Officers","Rape Report Rate","Sheriff",
                                                                    "Community Policing Component","Violent Crime Rate","Property Crime Rate"),
          label="femalecopsallarrestrate",
          title="The Effect of Female Police on All Crimes Arrest Rates, 1987-2016",
          dep.var.labels="All Crimes Arrest Rate",
          style="apsr",notes="Year dummy variables included.",out="./Female Cops Paper/Tables/femalecopsallarrestrate.tex",
          font.size = "small",table.placement="h!")

#table 5
model1 <- felm(sumpropcrimesrate~percfemale
               |year|0|0,data=genderpolice)
model1controls <- felm(sumpropcrimesrate~percfemale+totalemployees+log(budget)+percblackfulltime+propertycrimerate+sheriff+communitypolicing
                       |year|0|0,data=genderpolice) 
model2controls <- felm(sumpropcrimesrate~percfemale+totalemployees+ percblackfulltime+propertycrimerate+sheriff
                       |year|0|0,data=genderpolice)
model1a <- felm(propertycrimerate~percfemale|year|0|0,data=genderpolice)
model1acontrols <- felm(propertycrimerate~percfemale+totalemployees+log(budget)+percblackfulltime+sheriff+communitypolicing
                        |year|0|0,data=genderpolice) 
model2acontrols <- felm(propertycrimerate~percfemale+totalemployees+percblackfulltime+sheriff|year|0|0,data=genderpolice)
stargazer(model1a,model1,model2acontrols,model2controls,model1acontrols,model1controls,covariate.labels = c("Perc. Female Officers",
                                                                                                            "Total Agency Employees","Logged Agency Budget",
                                                                                                            "Perc. Black Officers","Property Crime Report Rate","Sheriff",
                                                                                                            "Community Policing"),
          label="femalecopspropcrimerate",
          title="The Effect of Female Police on Property Crime Arrest and Report Rates, 1987-2016",
          dep.var.labels=c("PC Report Rate","PC Arrest Rate","PC Report Rate",
                           "PC Arrest Rate","PC Report Rate","PC Arrest Rate"),
          style="apsr",notes="Year dummy variables included.",
          out= "./Female Cops Paper/Tables/femalecopspropcrimerate.tex",
          font.size = "tiny",table.placement="h!",column.sep.width = "-10pt")

#table 6
model1 <- felm(sumviolarrestsnoraperate~percfemale
               |year|0|0,data=genderpolice)
model1controls <- felm(sumviolarrestsnoraperate~percfemale+totalemployees+log(budget)+percblackfulltime+violentcrimerate+sheriff+communitypolicing
                       |year|0|0,data=genderpolice) 
model2controls <- felm(sumviolarrestsnoraperate~percfemale+totalemployees+ percblackfulltime+violentcrimerate+sheriff
                       |year|0|0,data=genderpolice)
model1a <- felm(violentcrimeratenorape~percfemale|year|0|0,data=genderpolice)
model1acontrols <- felm(violentcrimeratenorape~percfemale+totalemployees+log(budget)+percblackfulltime+sheriff+communitypolicing
                        |year|0|0,data=genderpolice) 
model2acontrols <- felm(violentcrimeratenorape~percfemale+totalemployees+percblackfulltime+sheriff|year|0|0,data=genderpolice)
stargazer(model1a,model1,model2acontrols,model2controls,model1acontrols,model1controls,covariate.labels = c("Perc. Female Officers",
                                                                                                            "Total Agency Employees","Logged Agency Budget",
                                                                                                            "Perc. Black Officers","Violent Crime Report Rate","Sheriff",
                                                                                                            "Community Policing"),
          label="femalecopsviolentcrimeratenorape",
          title="The Effect of Female Police on Violent Crime Arrest and Report Rates (No Rape Included), 1987-2016",
          dep.var.labels=c("VC Report Rate","VC Arrest Rate","VC Report Rate",
                           "VC Arrest Rate","VC Report Rate","VC Arrest Rate"),
          style="apsr",notes="Year dummy variables included.",
          out= "./Female Cops Paper/Tables/femalecopsviolentcrimeratenorape.tex",
          font.size = "tiny",table.placement="h!",column.sep.width = "-10pt") 

#table 7
main_model(genderpolice,"The Effect of Female Police on Rape Report and Clearance Rates, 1987-2016",
           "femalecopsmaincaseclosure",genderpolice$rapeclearancerate,genderpolice$raperate,genderpolice$percfemale,
           "Perc. Female Officers",
           c("Rape Report Rate","Rape Clearance Rate","Rape Report Rate",
             "Rape Clearance Rate","Rape Report Rate","Rape Clearance Rate"))

#table 8
model1 <- felm(sumrapearrestsrate~percfemale
               |STATECODE+year|0|STATECODE,data=genderpolice)
model1controls <- felm(sumrapearrestsrate~percfemale+
                         totalemployees+log(budget)+
                         percblackfulltime+raperate+sheriff+communitypolicing
                       |STATECODE+year|0|STATECODE,data=genderpolice)
model2controls <- felm(sumrapearrestsrate~percfemale+
                         totalemployees+
                         percblackfulltime+raperate+sheriff
                       |STATECODE+year|0|STATECODE,data=genderpolice)
model1a <- felm(raperate~percfemale
                |STATECODE+year|0|STATECODE,data=genderpolice)
model1acontrols <- felm(raperate~percfemale+
                          totalemployees+log(budget)+
                          percblackfulltime+sheriff+communitypolicing
                        |STATECODE+year|0|STATECODE,data=genderpolice)
model2acontrols <- felm(raperate~percfemale+
                          totalemployees+
                          percblackfulltime+sheriff
                        |STATECODE+year|0|STATECODE,data=genderpolice)
stargazer(model1a,model1,model2acontrols,model2controls,model1acontrols,model1controls,covariate.labels = c("Perc. Female Officers",
                                                                                                            "Total Agency Employees","Logged Agency Budget",
                                                                                                            "Perc. Black Officers","Rape Report Rate","Sheriff",
                                                                                                            "Community Policing Component"),
          label="femalecopsstatefixedeffectsandstateclusters",
          title="The Effect of Female Police on Rape Report and Arrest Rates, 1987-2016: Adding State Fixed Effects and State Clustered Standard Errors",
          dep.var.labels=c("Rape Report Rate","Rape Arrest Rate","Rape Report Rate",
                           "Rape Arrest Rate","Rape Report Rate","Rape Arrest Rate"),
          style="apsr",notes="State and year dummy variables included. SEs clustered by state.",out="./Female Cops Paper/Tables/femalecopsstatefixedeffectsandstateclusters.tex",
          font.size = "tiny",table.placement="h!",column.sep.width = "-10pt")

#table 9
model1 <- felm(sumrapearrestsrate~percfemale
               |fips+year|0|fips,data=genderpolice)
model1controls <- felm(sumrapearrestsrate~percfemale+
                         totalemployees+log(budget)+
                         percblackfulltime+raperate+sheriff+communitypolicing
                       |fips+year|0|fips,data=genderpolice)
model2controls <- felm(sumrapearrestsrate~percfemale+
                         totalemployees+
                         percblackfulltime+raperate+sheriff
                       |fips+year|0|fips,data=genderpolice)
model1a <- felm(raperate~percfemale
                |fips+year|0|fips,data=genderpolice)
model1acontrols <- felm(raperate~percfemale+
                          totalemployees+log(budget)+
                          percblackfulltime+sheriff+communitypolicing
                        |fips+year|0|fips,data=genderpolice)
model2acontrols <- felm(raperate~percfemale+
                          totalemployees+
                          percblackfulltime+sheriff
                        |fips+year|0|fips,data=genderpolice)
stargazer(model1a,model1,model2acontrols,model2controls,model1acontrols,model1controls,covariate.labels = c("Perc. Female Officers",
                                                                                                            "Total Agency Employees","Logged Agency Budget",
                                                                                                            "Perc. Black Officers","Rape Report Rate","Sheriff",
                                                                                                            "Community Policing Component"),
          label="femalecopscountyfixedeffectsandclusters",
          title="The Effect of Female Police on Rape Report and Arrest Rates, 1987-2016: Adding County Fixed Effects and County Clustered Standard Errors",
          dep.var.labels=c("Rape Report Rate","Rape Arrest Rate","Rape Report Rate",
                           "Rape Arrest Rate","Rape Report Rate","Rape Arrest Rate"),
          style="apsr",notes="County and year dummy variables included. SEs clustered by county.",out="./Female Cops Paper/Tables/femalecopscountyfixedeffectsandclusters.tex",
          font.size = "tiny",table.placement="h!",column.sep.width = "-10pt")

#table 10
fe <- genderpolice %>% 
  group_by(AGENCY) %>%
  filter(n()>1)
model1 <- felm(sumrapearrestsrate~percfemale
               |AGENCY+year|0|AGENCY,data=fe)
model1controls <- felm(sumrapearrestsrate~percfemale+
                         totalemployees+log(budget)+
                         percblackfulltime+raperate+sheriff+communitypolicing
                       |AGENCY+year|0|AGENCY,data=fe)
model2controls <- felm(sumrapearrestsrate~percfemale+
                         totalemployees+
                         percblackfulltime+raperate+sheriff
                       |AGENCY+year|0|AGENCY,data=fe)
model1a <- felm(raperate~percfemale
                |AGENCY+year|0|AGENCY,data=fe)
model1acontrols <- felm(raperate~percfemale+
                          totalemployees+log(budget)+
                          percblackfulltime+sheriff+communitypolicing
                        |AGENCY+year|0|AGENCY,data=fe)
model2acontrols <- felm(raperate~percfemale+
                          totalemployees+
                          percblackfulltime+sheriff
                        |AGENCY+year|0|AGENCY,data=fe)
stargazer(model1a,model1,model2acontrols,model2controls,model1acontrols,model1controls,covariate.labels = c("Perc. Female Officers",
                                                                                                            "Total Agency Employees","Logged Agency Budget",
                                                                                                            "Perc. Black Officers","Rape Report Rate","Sheriff",
                                                                                                            "Community Policing Component"),
          label="femalecopsfixedeffects",
          title="TThe Effect of Female Police on Rape Report and Arrest Rates, 1987-2016: Fixed Effects",
          dep.var.labels=c("Rape Report Rate","Rape Arrest Rate","Rape Report Rate",
                           "Rape Arrest Rate","Rape Report Rate","Rape Arrest Rate"),
          style="apsr",notes="Agency and year dummy variables included. SEs clustered by agency.",out="./Female Cops Paper/Tables/femalecopsfixedeffects.tex",
          font.size = "tiny",table.placement="h!",column.sep.width = "-10pt")

#table 11
model1 <- felm(sumrapearrestsrate~percfemale+percfemalesq
               |year|0|0,data=genderpolice)
model1controls <- felm(sumrapearrestsrate~percfemale+percfemalesq+
                         totalemployees+log(budget)+
                         percblackfulltime+raperate+sheriff+communitypolicing
                       |year|0|0,data=genderpolice) 
model2controls <- felm(sumrapearrestsrate~percfemale+percfemalesq+
                         totalemployees+
                         percblackfulltime+raperate+sheriff
                       |year|0|0,data=genderpolice)
model1a <- felm(raperate~percfemale+percfemalesq
                |year|0|0,data=genderpolice)
model1acontrols <- felm(raperate~percfemale+percfemalesq+
                          totalemployees+log(budget)+
                          percblackfulltime+sheriff+communitypolicing
                        |year|0|0,data=genderpolice)
model2acontrols <- felm(raperate~percfemale+percfemalesq+
                          totalemployees+
                          percblackfulltime+sheriff
                        |year|0|0,data=genderpolice)
stargazer(model1a,model1,model2acontrols,model2controls,model1acontrols,model1controls,covariate.labels = c("Perc. Female Officers","Perc. Female Officers Sq.",
                                                                                                            "Total Agency Employees","Logged Agency Budget",
                                                                                                            "Perc. Black Officers","Rape Report Rate","Sheriff",
                                                                                                            "Community Policing Component"),
          label="femalecopsnointeractionsquared",
          title="The Effect of Female Police on Rape Report and Arrest Rates, 1987-2016: Adding a Squared Term",
          dep.var.labels=c("Rape Report Rate","Rape Arrest Rate","Rape Report Rate",
                           "Rape Arrest Rate","Rape Report Rate","Rape Arrest Rate"),
          style="apsr",notes="Year dummy variables included.",out="./Female Cops Paper/Tables/femalecopsnointeractionsquared.tex",
          font.size = "tiny",table.placement="h!",column.sep.width = "-10pt")

#table 12
not100 <- genderpolice %>% filter(percfemale<100)

main_model(not100,"The Effect of Female Police on Rape Report and Arrest Rates, 1987-2016: Excluding All-Female Agencies",
           "femalecopsnooutliers",not100$sumrapearrestsrate,not100$raperate,not100$percfemale,
           "Perc. Female Officers",
           c("Rape Report Rate","Rape Arrest Rate","Rape Report Rate",
             "Rape Arrest Rate","Rape Report Rate","Rape Arrest Rate"))

#table 13
dta_m <- match.data(mod_match) 
model1_match <- felm(sumrapearrestsrate~percfemale
                     |year|0|0,data=dta_m) 
model2controls_match <- felm(sumrapearrestsrate~percfemale+
                               totalemployees+
                               percblackfulltime+raperate+sheriff
                             |year|0|0,data=dta_m)
model1a_match <- felm(raperate~percfemale
                      |year|0|0,data=dta_m) 
model2acontrols_match <- felm(raperate~percfemale+
                                totalemployees+
                                percblackfulltime+sheriff
                              |year|0|0,data=dta_m)
stargazer(model1a_match,model1_match,model2acontrols_match,model2controls_match,covariate.labels = c("Female Officers Over 8 Perc.",
                                                                                                     "Total Agency Employees",
                                                                                                     "Perc. Black Officers","Rape Report Rate","Sheriff"),
          label="femalecopsmatching",
          title="TThe Effect of Female Police on Rape Report and Arrest Rates, 1987-2016: Matching Analysis",
          dep.var.labels=c("Rape Report Rate","Rape Arrest Rate","Rape Report Rate","Rape Arrest Rate"),
          style="apsr",notes="Year dummy variables included.",out="./Female Cops Paper/Tables/femalecopsmatching.tex",
          font.size = "tiny",table.placement="h!",column.sep.width = "-10pt")

#table 14
pre1993 <- genderpolice %>% filter(year<1993) 
model1 <- felm(sumrapearrestsrate~percfemale
               |year|0|0,data=pre1993)
model1controls <- felm(sumrapearrestsrate~percfemale+totalemployees+
                         percblackfulltime+raperate+sheriff
                       |year|0|0,data=pre1993)
model2controls <- felm(sumrapearrestsrate~percfemale+
                         totalemployees+
                         percblackfulltime+raperate+sheriff
                       |year|0|0,data=pre1993)
model1a <- felm(raperate~percfemale
                |year|0|0,data=pre1993)
model1acontrols <- felm(raperate~percfemale+
                          totalemployees+
                          percblackfulltime+sheriff
                        |year|0|0,data=pre1993)
model2acontrols <- felm(raperate~percfemale+
                          totalemployees+
                          percblackfulltime+sheriff
                        |year|0|0,data=pre1993)
stargazer(model1a,model1,model2acontrols,model2controls,model1acontrols,model1controls,covariate.labels = c("Perc. Female Officers",
                                                                                                            "Total Agency Employees",
                                                                                                            "Perc. Black Officers","Rape Report Rate","Sheriff"),
          label="femalecops1997before",
          title="The Effect of Female Police on Rape Report and Arrest Rates, 1987-1993",
          dep.var.labels=c("Rape Report Rate","Rape Arrest Rate","Rape Report Rate",
                           "Rape Arrest Rate","Rape Report Rate","Rape Arrest Rate"),
          style="apsr",notes="Year dummy variables included.",out="./Female Cops Paper/Tables/femalecops1997before.tex",
          font.size = "tiny",table.placement="h!",column.sep.width = "-10pt")

#table 15
post1993 <- genderpolice %>% filter(year>1993)

main_model(post1993,"The Effect of Female Police on Rape Report and Arrest Rates, 1997-2016",
           "femalecops1997on",post1993$sumrapearrestsrate,post1993$raperate,post1993$percfemale,
           "Perc. Female Officers",
           c("Rape Report Rate","Rape Arrest Rate","Rape Report Rate",
             "Rape Arrest Rate","Rape Report Rate","Rape Arrest Rate"))

#table 16
gp_1987 <- genderpolice %>% filter(year==1987)
gp_1990 <- genderpolice %>% filter(year==1990)
gp_1993 <- genderpolice %>% filter(year==1993)
gp_1997 <- genderpolice %>% filter(year==1997)
gp_2000 <- genderpolice %>% filter(year==2000)
gp_2003 <- genderpolice %>% filter(year==2003)
gp_2007 <- genderpolice %>% filter(year==2007)
gp_2013 <- genderpolice %>% filter(year==2013)
gp_2016 <- genderpolice %>% filter(year==2016)

stargazer(year_fun_arrests_pre(gp_1987),year_fun_arrests_pre(gp_1990),year_fun_arrests_pre(gp_1993),
          year_fun_arrests_pre(gp_1997),year_fun_arrests_post(gp_2000),year_fun_arrests_post(gp_2003),
          year_fun_arrests_post(gp_2007),year_fun_arrests_post2007(gp_2013),year_fun_arrests_post2007(gp_2016),covariate.labels = c("Perc. Female Officers",
                                                                                                                                    "Total Agency Employees","Logged Agency Budget","Perc. Black Officers","Rape Report Rate","Sheriff","Community Policing","Domestic Violence Unit"),
          label="femalecopsarrestsbyyear",
          title="The Effect of Female Police on Rape Arrest Rates, 1987-2016: Regressions By Year",
          column.labels=c("1987","1990","1993","1997","2000","2003","2007","2013","2016"),
          style="apsr",out="./Female Cops Paper/Tables/femalecopsarrestsbyyear.tex",omit.stat=c("ser"),
          font.size = "tiny",table.placement="h!",dep.var.labels = "Rape Arrest Rate")

#table 17
stargazer(year_fun_reports_pre(gp_1987),year_fun_reports_pre(gp_1990),year_fun_reports_pre(gp_1993),
          year_fun_reports_pre(gp_1997),year_fun_reports_post(gp_2000),year_fun_reports_post(gp_2003),
          year_fun_reports_post(gp_2007),year_fun_reports_post2007(gp_2013),year_fun_reports_post2007(gp_2016),covariate.labels = c("Perc. Female Officers",
                                                                                                                                    "Total Agency Employees","Logged Agency Budget","Perc. Black Officers","Sheriff","Community Policing","Domestic Violence Unit"),
          label="femalecopsreportsbyyear",
          title="The Effect of Female Police on Rape Report Rates, 1987-2016: Regressions By Year",
          column.labels=c("1987","1990","1993","1997","2000","2003","2007","2013","2016"),
          style="apsr",out="./Female Cops Paper/Tables/femalecopsreportsbyyear.tex",omit.stat=c("ser"),
          font.size = "tiny",table.placement="h!",column.sep.width = "-5pt",dep.var.labels = "Rape Report Rate")

#table 18
sheriff <- genderpolice %>% filter(sheriff==1)
main_model(sheriff,"The Effect of Female Police on Rape Report and Arrest Rates, 1987-2016: Only Sheriffs' Offices",
           "sheriffspecification",sheriff$sumrapearrestsrate,sheriff$raperate,sheriff$percfemale,
           "Perc. Female Officers",
           c("Rape Report Rate","Rape Arrest Rate","Rape Report Rate",
             "Rape Arrest Rate","Rape Report Rate","Rape Arrest Rate"))

#table 19
municipal <- genderpolice %>% filter(municipal==1)
main_model(municipal,"The Effect of Female Police on Rape Report and Arrest Rates, 1987-2016: Only Municipal Police Departments",
           "municipalspecification",municipal$sumrapearrestsrate,municipal$raperate,municipal$percfemale,
           "Perc. Female Officers",
           c("Rape Report Rate","Rape Arrest Rate","Rape Report Rate",
             "Rape Arrest Rate","Rape Report Rate","Rape Arrest Rate"))

#table 20
blackandcity <- genderpolice %>% filter(year==2016 | year==2013 | year==2007 | year==2000 | year==1990)

model1 <- felm(sumrapearrestsrate~percfemale
               |year|0|0,data=blackandcity)
model1controls <- felm(sumrapearrestsrate~percfemale+
                         totalemployees+log(budget)+
                         percblack+percblackfulltime+log(1+population)+raperate+sheriff+communitypolicing
                       |year|0|0,data=blackandcity)
model2controls <- felm(sumrapearrestsrate~percfemale+
                         totalemployees+log(budget)+
                         percblack+percblackfulltime+raperate+sheriff
                       |year|0|0,data=blackandcity)
model1a <- felm(raperate~percfemale
                |year|0|0,data=blackandcity)
model1acontrols <- felm(raperate~percfemale+
                          totalemployees+log(budget)+
                          percblack+percblackfulltime+sheriff+communitypolicing
                        |year|0|0,data=blackandcity)
model2acontrols <- felm(raperate~percfemale+
                          totalemployees+log(budget)+
                          percblack+percblackfulltime+sheriff
                        |year|0|0,data=blackandcity)
stargazer(model1a,model1,model2acontrols,model2controls,model1acontrols,model1controls,covariate.labels = c("Perc. Female Officers",
                                                                                                            "Total Agency Employees","Logged Agency Budget",
                                                                                                            "Perc. Black","Perc. Black Officers","Logged Population","Rape Report Rate","Sheriff",
                                                                                                            "Community Policing Component"),
          label="femalecopsaddcontrols",
          title="The Effect of Female Police on Rape Report and Arrest Rates, 1987-2016: Adding Place Controls",
          dep.var.labels=c("Rape Report Rate","Rape Arrest Rate","Rape Report Rate",
                           "Rape Arrest Rate","Rape Report Rate","Rape Arrest Rate"),
          style="apsr",notes="Year dummy variables included.",out="./Female Cops Paper/Tables/femalecopsaddcontrols.tex",
          font.size = "tiny",table.placement="h!",column.sep.width = "-10pt")

